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TV destinations in the plane {Pj : j = 1, . . ., N} are given as independent random variables with 
specified probability densities, and the problem is to find the location of the point P which minimizes 
the expected sum of the Euclidean distances PPj, 

In this paper, upper bounds for the minimizing sum of distances are found in terms of solutions to 
corresponding deterministic problems and the first and second moments of the probability densities. 
Three commonly occurring classes of bivariate probability densities: (A) normal. (B) exponential, and 
(C) symmetric exponential, are then considered. Numerical tests are presented which show that in all 
cases, Steffensen's iteration is effective in accelerating convergence. Finally it is shown that in contrast 
to the deterministic case, P need not be in the convex hull of the means of Pj and a sufficient condition 
is given for P to be in this convex hull. 

Key words: Always convergent algorithms; exponentially distributed points: facility location; location 
theory: normally distributed points: numerical analysis: optimization; stochastic Weber problem. 

1. Introduction 

A probabilistic extension of the classical Weber (Steiner) problem was introduced in [3|' and 
[6|. /V destinations in the plane, {Pj :j=l, . . ., N} , are given as random variables witb specified 
probability density functions. It is desired to find the location of a point P with ("artesian coordinates 
a= (a, ft) which minimizes the expected sum of the Euclidean distances PPj. A possible application 
of this problem to minimizing capital costs in water quality management is given in [6]. For sim- 
plicity we restrict ourselves in what follows toE 2 . 

Specifically, let P, be a random variable with probability density function </>j(#j), xj= (xj, 
7j),7=l,. . . , /V, where Xj e Sj C E 1 . Sj is the sample space of Pj. Denote by |jc — jcj| = [(% — Xj) 2J r 
(y — jj) 2 ] 1/2 the Euclidean distance between x= (x, y) and Xj. The expected value of PPj is then 

E ( \x — Xj | ) = I I \x — Xj | 4>j (xj ) dxj dyj and we wish to minimize the objective function 

4j(x) = S\ (3j \ \x-Xj\cpj(xj)dxjdyj (1) 

over Jc, where theft are given, positive weights. 

In this paper Sj = S for 7=1, . . . , N; that is, all Pj are in the same sample space. We also 
restrict P to S. It was shown in [6| under mild assumptions on <f>j that ifj{x) is strictly convex, and 
that the minimizing point P is unique. An iterative algorithm for finding P was obtained formally 
by setting Vijj(a) =0, by "solving" this equation in the form a = //(a), and then by passing to the 
iteration scheme x n+l = H(x n ) which this suggests. Specifically, 
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where £ = (£, 17) and G(£) = V j3/0j (f). The algorithm was shown in [6] to be a descent method 

j 
which is globally convergent to P and convergence was shown to be linear, locally. Although the 
considerations in [6] are for S = E 2 , the same proofs apply with minor variations when S= [0, 00) 
x [0, 00). 

In this paper, we first obtain upper bounds for 



i// = min \fj(x) 

xeS 



(3) 



in terms of the solutions to corresponding deterministic problems with the new destinations at the 
means of the random variables Pj. It is of interest to note that these bounds involve only the first 
and second moments of (f)j(xj). We also give a sufficient condition for P to be contained in an as- 
sociated convex hull. This is important because whereas in the classical deterministic problem P 
is always in the convex hull of the destinations (see [5] and [7]), a corresponding statement for the 
probabilistic problem is not necessarily true, as we demonstrate later with a counterexample. 
Next we consider three specific classes of probability density functions: 



(A) S = £2 and (/),(*., 



\27TO-j 



7- Vi 



# 



exp 



2VT 



Pi 
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■2 A , fJZZa) (*!Z»u) + (vzaii 



the bivariate normal density function, 
(B) S=[0, oo)x[0, °°),and 

</>./ (xj ) = kxjkyj exp {-XxjXj - k yj yj } , 
the bivariate exponential density function, and 
(C) S = £ 2 and 

</>j (xj ) = ~ k Xj k yj e xp {- X Xj \xj - Xj |- X yj I yj - Vj | } , 

a bivariate exponential density function with mean at {d x ., 6 y .) and standard deviation (V2/X^., 
y/2lky. ). In each case we seek the minimizing xeS. 

A major part of this paper deals with a special case of (A), namely the symmetric bivariate 
normal density when Xj and yj are independent (pj=0), and <jy.— (Ty* = ctj, /= 1., . . . , N. In 
this case, which is of considerable practical importance, all integrals appearing in (1) and (2) 
can be evaluated in closed form in terms of products of exponentials and modified Bessel Functions 
of orders and 1. This may be very useful in practical computations, since tables or subroutines can 
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then he used instead of time-consuming numerical quadrature. We also prove that in this case P is 
in the convex hull of the means of the destinations /jl,= (jjlx , l^u ) » which we denote by co{/x;}. 
By means of a simple counterexample with N=2 we show that even if p y = 0, when <j x .^cry. 
I* need not be in co {fij}. In the general case in (A) where the density is asymmetric, all double 
integrals are reducible to single integrals. This is also true for (B) and (C). 

Finally, it was suggested in [6], that since convergence is locally linear in (2), Steffensen's 
iteration [4] for acceleration of convergence should be effective. We report the results of randomly 
generated numerical examples which confirm that this is indeed the case. We also compare our 
upper bound for i// with its true value for some numerical examples. 



2. Upper Bounds for i/j 

An upper bound for (3) follows from Schwarz's inequality. Since <f)j(xj) ^ 0, and 

I I (f>j(xj)dxjdyj= 1 , 

we have 

i//U, y) = ^ ft f f VU-f) 2 + (y-v) 2 <Mf , y)d£dri 
j J J s 

ft <Mf, r))dedv] 112 



2#/{ f [ [U-0 2 +(y-^) 2 ]<M^)^7? 

j L J J s 
£ft x*-2xfi x . f£[*J]+y 2 -2y^ y . + £[yJ] 



S& 



( x — fi. rj ) 2 + ( v* — M yj ) 2 + o- 2 . . + o~l . 



1/2 



where /u^.= i?[jtj] , /x,y.== Z?[y/] are the respective jc and y means of c/> 7 - (jc; ) , and where 
a' 1 = £[%?] - M 2 . , (7 2 = ETyf] — fli are the respective x and y variances of 4>j(xj). This gives 

],^mmYflj[(x- f i x .) 2 + (y- /*,,)* + a*] 1 /* (4) 

where cr 2 = crl + o\ 2 ; . A coarser but more easily calculated upper bound follows immediately 

J x j y J 

from (4): 

^ min V Pj ([(x- il x ) 2 + (y- /x y . ) 2 ]'/ 2 + cr,) 

j 

2 jSjffj I ...in ^ Pj [ U - M.- ; ) 2 + (y- M,;) 2 ]" 2 = * Bound. (5) 

j j 

The second term on the right in (5) is the minimum weighted sum of the distances in a deterministic 
problem with the destinations at the means of Pj , and the first term is a weighted sum of the stand- 
ard deviations of the Pj. Note that for any 4>j(xj) n the bounds in (4) and (5) depend only upon its 
first and second moments. 
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3. The Bivariate Normal Density Function 

Assume that (f>j(xj) for 7= 1, . . ., N is given by (A). We now show that 

if pj = and o\ r . = cr y . = a , for 7=1, . . ., N then ~ae co {jx, } N = i (£\ 



In this case 



<M*))= (27ro-. 2 )- 1 exp{-^ [(*,-^.) 2 + (y.;-/x,,.) 2 ]}. (7) 

We first prove two lemmas which are useful in locating a. 

Lemma 1: Suppose z x (y) are differentiable, real valued functions defined on E 2 , i= 1, . . ., N and 
/e£ y be a point at which z(y) = X Zj(y) is minimized. Then y e co {yil^Lj 15 assured under the following 
two further assumptions: 

(a) The level curves S/={y e E 2 |zi(y)= c} are concentric circles centered at yj, £/ia£ is,Sj = 

{y||y-yil = r(c)}. 

(b) For some k, Vz k (y) # i/y ^ y k . 

PROOF: Assume, to the contrary, that y ^ co{yj}f=i. Then there is a line / containing 7, with a 
normal a* such that all fu i — 1 , . . . , iV are contained in an open half plane bounded by /. For 
y on /, we have a -y = a -y = 6, and a ■ (ft — y) > 0, for £ = 1, . . . , N. Now, because of (a), 

V*0)0 = |Vzi(y)| J 



l7-y/| 

and because of (b) | Vza (y ) | > 0. Therefore, we have 

(Vz(y))-a = y |Vz,(?) | j^- -a < 

i It — r/| 

which implies that z(y) is decreasing in the direction a from y. This contradicts the assumption 
that z(f) is minimized at y, hence y € co{y;}^ =1> 

LEMMA 2: Suppose that the probability density function (f)^(xj) is radially symmetric about its 
mean jjl } ; that is, suppose (/)j (xj) = f j ( | x ; — /x j | ) , x- } e E 2 . 7 T &e/i 

zj(i3i=i8 J |J E J*-iU J (?)dfdi, (8) 

satisfies assumption (a) in lemma 7, with y j = /Xj. 

PROOF: Use polar coordinates £ — /x x . = p cos 0, 77 — fiy.= p sin in (8) to obtain 

zj(x) = iijjj^ \x-i \<i>j(i)dZdT, = pjf f^ \x-i\M\£ r iLj\)d€dri 

= 1 I p [(* — /t*, — p cos 0) 2 + (y-ixy.-p sin6) 2 ] " 2 fj(p)dOdp 

= j o " pfj(p)dp j*' [r 2 + p 2 - 2pr, cos (0 - 4>j)]^ dQ (9) 

where r? = (x — p ^ ) 2 + (y — p tf ) 2 , 
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cos<|>,= (x-ii Xj )lr j9 sin<&j= (y - /ul !I} ) lr h (10) 

Using the periodicity of the integrand over [0, 2rr] , and letting y \ r = — <$>,, (9) becomes 

z,(*) = f*pfj(p)dp f^ [rf + p 2 -2prj cos ^^dV (11) 

Jo Jo J 

which implies that zj(x) in its dependence upon x is a function of rj, The Euclidean distance to 
/Xj, alone. Hence, the level curves of z/(#) are concentric circles centered about jlj. This proves 
the lemma. 

Now we wish to apply lemma 1 with 



Zj(*) = /BjJ"Jj*-?|to 



(£)d£dT) 



and with 0, (£) given by (7). It follows from [5] with /V = 1 that zj (x) is strictly convex on E 2 . Also 

V 2 , (x ) = ft j j ||^|j </>,• (| ) <#*,,, 

and an easy computation shows that Vzj(jZj) = where jLj = (M/., Pu)- Since Z/(Je) is strictly 
convex, Vz,- (* ) ^ if x ¥" JIj i for allj= 1, . . .,/V. Hence Zj(Jc) satisfies assumption (b) of lemma 1. 
In order to show that zj{x) satisfies assumption (a) of lemma I, observe tbat from (7) we have 



if)j (xj) = (2mrj ) ■ exp j -^—7 \xj - fij 



so that (f)j{xj) — fj (\x; ~ /I/|). By lemma 2, then, 2; (x) satisfies condition (a) of lemma 1, with 
Yj = juij. Now apply lemma 1 with y;= /Z ; - to obtain (6). 

The result in (6) is not trivial. In the classical deterministic problem 5 is always contained in 
co{Pj}y =l where Pj are the given destinations (see [5] and [7] ). However, in the probabilistic 
case, we now show by means of a simple counterexample with N = 2, that even when the </>; are 
normal bivariate density functions with pj — 0, if <j x . ¥^ <j u . then a is not necessarily in co{/Z/ }f ssV 

Consider the problem: N = 2, (£/(#/) for j = 1, 2 given by (A), 

/3i = ^ 2 = 1, Mxi = 1? /^y^ 0' 0"*i = cr, (r yi = 1; /Aa- 2 = 0>/Xy 2 = 1, cr*«== 1, CTy 2 = o-: 

pj = 0. In this problem co{/Z ; } 2 = 1 — { (x, y) \x + y = 1,0^x^1}, the closed line segment joining 
the means. From symmetry considerations, it is easy to show that for each fixed cr, < cr < 00, 

the minimizing point a = (a, (3 ) satisfies a — (3. Therefore a is in co{/Zj }r ^ if and only if a = (3 = -. 

It also follows using Vi//(«) =0 that a satisfies — (a, a) = for each cr: 

ax 

+ /J..-=3-'H(^ + '' , M-* 
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We give the following argument: As a -> it is well known (see [2], p. 324 and [8], chapter 1, 
for example) that in the sense of distributions (V 2tt a) l exp {— 1/2 (u— l)' 2 /cr 2 } -> 8(u — 1), the 
Dirac-delta function centered at u= 1. In each integral of (12), one dimension of the probabilistic 
problem then becomes deterministic with the corresponding coordinate of the destination fixed 
at 1. Lettingor -» in (12) we obtain 



n<X} J. V(a-f) 2 +(a-D 2 ^ J.V(«-1)H(«-^ 



= I' 2a " f " 1 ■ e-«V^0. (13) 

J V( a -^) 2 +(«-l) 2 ? 



We now show that the unique a satisfying (13), also satisfies 1/2 < a < 1. Continuity arguments then 
imply that for cr sufficiently small the unique a satisfying (12) is also such that 1/2 < a < 1, i.e., for 
cr sufficiently small a 4 co{ju;}j =l , thereby providing a counterexample. We have 



'©"J" , , '., '-"* 



"iN 



/: 




df<0 



and 



^ (1)= / 3 ljT^ e ^ f ^ = /_ 1 1 ei '^ >0 - 



Hence a satisfying (12) is such that ^ < a < 1 , which concludes the proof. (In fact a —* 0.66 approxi- 
mately, as cr^ 0, as seen from table 1.) 

For purposes of illustration we have solved (12) numerically using an iterative scheme analogous 
to (2). In table 1 we give a for various values of cr, < cr < oo. As cr goes from to oc , a goes mono- 
tonically from 0.660 to 0. Therefore for no value of cr except cr= 1 (when cr Xl = cr Ul = 1 = cr.r 2 = <j lhz ) 
is a in co{/X y }] =1 . This is in marked contrast to the deterministic case when a is always in the 
convex hull of the destinations. We also give in table 1, i//, the bound for \fj from equations (4), 
\\t BOUND from (5), and M, the number of iterations of (2) to achieve three place accuracy. The 
bound for i// from eq (4) was calculated from an obvious extension of the classical iteration for solv- 
ing the Weber problem (see [5], or [6], for example). 
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TABLE 1. Counterexample 



IT 


a 


M 


«/' 


Bound for.// 
from (4) 


./; hound 

from (5) 


1/32 


0.660 


7 


2.123 


2.450 


3.415 


1/16 


.658 


7 


2.125 


2.452 


3.418 


w 


.650 


7 


2.134 


2.462 


3.430 


L/4 


.623 


7 


2.173 


2.500 


3.476 


M2 


.567 


8 


2.326 


2.646 


3.650 


L* 


.500* 


— 


2.814* 


— 


4.243 


2 


.416 


9 


4.074 


4.690 


4.886 


4 


.391 


11 


6.969 


8.367 


9.660 


8 


.259 


11 


13.108 


16.186 


17.539 


16 


.122 


12 


25.708 


32.094 


33.477 


32 


.043 


12 


51.138 


64.047 


65.445 


64 


.012 


12 


102.136 


128.023 


129.430 


128 


.003 


12 


204.202 


256.012 


257.422 


256 


.001 


12 


408.370 


512.006 


513.418 


512 


.000 


12 


816.723 


1024.003 


1025.42 



^Computed from clos ed form expressions (see (24)). 
<£BOUND = V2 + 2Vl + o---. 

The minimum i//(x, y) = »/i mi „ in the corresponding deterministic 
problem with /', at jlj is: t// mill = Vz. 



3.1. Evaluation of Integrals in the Symmetric Uncorrelated Case 



We now give explicit closed form expressions for the evaluate 



of the Integrals in (2) 



hen p, = and cr., ■= cr y . = tr, for j= L, . . ., N. In this case <f>j(xj) is given by (7) and the right- 



hand side of (2) becomes 



H{x) 



2 &M(*) 

2 fijDjix) 



(14) 



NAx 



*) = (2iro-})-> r J |[U-|) 2 +(y-T,) 2 ]-" 2 exp 



2(7j 



(15) 



/) j (x) = (27ro-p- 1 | OC | X [U-^) 2 +(y~7 ? ) 2 ]- 1 / 2 exp{-^[(^- Mx p2 +( ^_ My p 2] J^ d ^ 

(16) 

We evaluate the first component of Nj(x), which for simplicity we denote by Nj(x). The evaluation 
of the second component follows similarly. Substitute £ — x = p cos 0, r\ — y = p sin to obtain 

= (2tto-j) ' I (jt + p cos 6) exp j --—[(x-u x . + p cos 6>) 2 

Jo Jo I 2o"7 ' y 



+ (y-/x /y . + p sinfl) 2 ] \d0dp = K j (x) J rL l (x). 
Use the notation in (10) and let 
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(17) 



Jj(p) = exp {-(prj/o-?) cos (0-Qj)}dO 
J o 

= I exp {- (prjlcrj) cos ^JcW 

where we have substituted — <&j= x P and used the periodicity of the integrand. Then, we have 

Kj(x)=x(2tt(t])^ exp {-rj/2crj} Pexp {- (p 2 l2<r1)}Jj(p)dp. 

Jo 

Reference to [1], p. 376, 9.6.16 or 9.6.19 shows that 

J j (p)=27rh(pr j l<TJ) 

where /& is the modified Bessel function of the first kind of order k. Reference to [1], p. 487, 11.4.31 
then gives 



Kj(x)= (xl*,)W2)* exp {-rj/4o-]}/o{r]/4o-?}. 
In order to evaluate Lj (x) note that 

Lj(x)= (2mrJ)-' exp {-r]l2<rj} \p exp {-pV2<7J}Mj(p)dp 

JO 

where 

f 27r 
Mj(p)= cos <9exp {(—prj/dj) cos (0-O;)}d0. 
Jo 

Let # — ^j— ^, use the peridoicity of the integrand and refer to [1], p. 376,9.6.19, to obtain 

M j (p)=-27T((x- f x Xj )lr j )I 1 {pr j /a]}. 
Substituting in (19), this gives 

Lj(p) = -((x-tL Xj )lrj) exp {-rj/2crj} o"/ 2 J p exp {- p 2 /2aj}I l {pr j laj}dp . 

Integrating by parts in the last integral, we obtain 

Lj(p)=-((x-fx Xj )/rj) exp {-rjl2crj}. 



(18) 



(19) 



■Idpnlor]} exp -p 2 /2o-j 



"+(rj/cr|) (""exp {-p 2 /2o-?}/;{po/cr|}rfpl. 

=o Jo J 



(20) 



The integrated term vanishes at p = since 7i(0)=0. At p=oo we use I\{z) — e z (2irz)~ 112 (see 
[1], p. 377, 9.7.1) to show that the integrated term vanishes there too. The expression in (20) can 
be further simplified by using 21 [ = I + 1> (see [1], p. 376, 9.6.26). This simplification combined with 
use of [1], p. 487, 11.4.31 gives 
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/., / (p) = -((^-^, j )/2o- j )( 7 7/2)>/ 2 exp {-/-5/4crJ}(/o{r5/4cr5}+/,{r5/4«r5}). (21) 

Substituting (18) and (21) in (17) gives the explicit closed form expression for a typical double 
integral in (14), 

Nj(X)=(2aj) l (*72) 1/2 exp {-r?/4c7?}(U + ^.)/o{r)/4o-]} - {x- n x .)h{r)lkr]))}. (22) 

Calculations completely analogous to those which gave Kj(x) in (18) also yield 

Dj(x) = oy 1 0r/2) 1 ' 1 exp {-r]/4o-?)/o{r]/4o-?}. (23) 

(22) and (23) are the closed form expressions required for the iterative algorithm (2), in this special 
case. 

The equation for the minimizing x, x = H(x) can now he written as 



2 2 PPf 1 ex P {-rjl^TJjiix^JZ^hirjI^rJj-ix-fl^hirjI^T]^ 

. 1 



(24) 



SiS/r^exp {-r?/4o-|}/o{rf/4o-?} 
j 

It is of interest to note that (24) can be rewritten in the equivalent form 

.v 
x = 2 Aa-m/,, (25) 



where 

£,(r,' exp {-r|/4o-|}(/ {r|/4o-|.}+/ 1 {r A 2 /4o- A 2 }) 



X/,=- y 

2 oyi exp {-r?/4o-?}(/ {r?/4o-?} + /,{r?/4o-?}) 

satisfy < X* < 1, EXa = 1. Therefore (25) is a closed form expression (in terms of x) for Jc as a 
convex combination of jut/, , with constants X/,. We proved earlier that in the case of bivariate normal 
density functions with Xj and yj independent and a x j= cr y ■= (Tj , such a combination always exists 
since xeco { ^} N k=zt . 

It is also possible to evaluate the objective function ty(x) in (1) without numerical quadrature. 
We have 



and 



Zj(*)=^(27ro-)) ' exp{-r?/2o-j}| pVj(p) exp {-p 2 /2o-?} rfp 

= /3.;o-j 2 exp {-r?/2o-?} | " p*/ (po/o-}) exp {-p2/2o- 2 ,.} dp. (26) 

JO 
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Refer to [1] p. 487, 11.4.31 and differentiate with respect to the parameter a 2 appearing there. 
This gives a closed form expression for the integral in (26), and (26) then becomes 

Zj(x)=((3 J l<T j )(nl2yH exp {-rf/4<rj} ((r|/2)(/ (r|/4o-|) 

+ /i(r|/4o-|))+(r|/o(r?/4o-p). 

3.2. Evaluation of Integrals in the General Case 

To conclude our analysis of the bivariate normal density function we now show that even in 
the general case (when the density is not symmetric and uncorrected), all of the double integrals 
appearing in (1) and (2) can be evaluated as single integrals in terms of error functions. The iteration 
in (2) then becomes considerably more efficient because only single integrals need be evaluated 
by numerical quadrature, using subroutines available for the calculation of the error function. 

In the general case (f)j(xj) is given by (A) and the right-hand side of (2) becomes 



H(x) 






(27) 



/" J"* 

J — x J — oc 



(%-£) 2 + (y-n) 2 ]" 1/2 exp 



Nf(x)=(2ir<r x fTyj VT^j 

1 



2vT 



"2pj 



Pj 



f-M*i 



CTj 



(T XJ 



(Tu 



y-v>uj 



(Ti 



d^dr] (28) 



D*(x)=(27rcr x /r yj VT^)" 1 . 

| X | X [U-^) 2 +(y~7|) 2 ]- 1 / 2 exp 



2vT 



-Pj 



£-***, 



CTj 



-2 P ,'^ 






d£dr). (29) 



Again, we evaluate only the first component of N* (x) which, for simplicity, we denote by /V* (x). 
The evaluation of the second component follows similarly. Substitute, as before, £ — x = p cos 0, 
tj — y= p sin 6. The exponent in brackets in (28) and (29) then becomes 



: — /x, r . + p cos 0\ 2 /x — pj .-\- p cos 0\ / y— /x y + p sin \ /y— /x y . + p sin 0\ 2 



**J 



(Tu 



R \ 2 /?- 
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wh( 






(x — fJLx) COS (x— u r ) sin 6+ (y—u,,.) cos (9 (j~~"My-) sni 8 
Bj = Bj(6) = -f- -f>> 2 ~ + U ^~ 



r 2 



x — fl., 






Note that /4 j = /4j(0) is never zero for < o>., <r ,, . < oo. /V*(*) now becomes 



/V*(*) = (27ro- J ./j-„ ,VT^5) ' exp {-rll(2VT=ri)} exp {fi«/(a4 J Vl^})}P J (fl)de 

J J In 



Pj(0)= (* + pcos 0) exp {-^j(p + Bj//l;) 2 /(2VF^)}dp. 

We now show thai Pj(6) can be evaluated in closed form. For each fixed #, lei 

BA/ A, \' 12 



t=(p+- 



Aj/VZVl-f 



lllrll 



P,(0) = (2VP^]//C)'/ 2 . 



(30) 



Bj 



---cos 6»1+ (2Vl- p ]Mj)" 2 (cos 0)t e (2 cfc 

J HI2.4.V \- p t) '» 

= (77Vr^i/Z4 j )'/ 2 (x-^cos 0) erfc {fij/^VF^j)'/*} 



+ (VT^M J )(cos 0) exp {-fl*/(2^VF^)} 

where erfc is the complementary error function (defined, for instance, in [1], p. 297). 
Similarly, for the terms in the denominator of (27), we have 



Df(jt)=(2™.,..<r, / .Vl-p;j) ■ exp {-r7/(2VT^pl)} exp {flJ/ZfjVT^WjW* 



(M0)= ["exp {-// J (p + S J /// J ) 2 /(2Vr=7i)K/p 

J<) 

= (ttVF^W 2 erfc {B j /(24 j Vf=^)^}. 
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Aslo, to evaluate i//(x) in (1) we have 



<K*)=2)zj(*) 



(31) 



t'here 



Zj(x)- 



■■f3A2ncr, r „ j VY^ l )-' j ^j J( x -^+ (y-v) 2 ] m {- 2 ^ x l _ ; 



^-MxA 2 



■2p,|^^V^^U^^^V 



0\r.. / \ o- 



1 ; /; 



o"«, 



dgdrj. 



This leads to 

T277- 

z J (*)=/8j(2mr» J «r,.vT := ^j)- 1 exp {-rJ7(2Vp=^5)} exp {B*,l(2AjVl=tf)}Rj(6)dd 

■' J ' Jo 

where 

Rj{0)=f*p> exp {-^j(p + BjMj)»/(2Vr=^J)}rfp. 

Use the substitution (30) to obtain 

R J (e) = (2VT^iA J )^r (- 

Integration by parts in the first term then gives 



— n 2 B ; /2Vl — n 2 \ l / 2 #?\ 



«j(fl) = (2Vl=pjMj)^ 



!*- PJ + ^~) erfc ^/(^j^v^^j)^} 



2^j V ^ , 



exp {- B»jl(2Aj VT=p])} 



4. The Bivariate Exponential Density Function 

In this case all double integrals appearing in (1) and (2) can be reduced to single integrals. 
Using the notation in (14) we have 

Nj(x) = \ Xj k y . f " f " f[ (* - f ) 2 + (y- r,) 2 ] -" 2 exp {-(A,,f + X yj i)) }<tfdij 

Dj(x)=k x .ky. j" J" [(*-£)* + (r-r ? ) 2 ]- 1 / 2 exp{-(X^ + X, / .T,)}^ < /7 ) . 

We consider only the first component of Nj(x) which we denote by Nj(x). Let £ — x = p cos 0, 

17 — y = p sin 0, r = x 2 -\-y 2 , cos = — #/f, sin = — y/f. Thus, tan = yjx and rr ^ ^ 3tt/2. 
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Then we obtain 
Nj(x) = \ X .X U . 



77-/2 f to /•() /•-.r/cosG T277- r- y/ sin e 

+ 1 I +11 ) ( U + pcosG) 



= ^A yj exp {— xk Xj — yk y 



r() r .r/oosB r 2tt r 
Jn/2 Jo Je Jo 

exp {— k Xj (x + pcos 0) — k yj (y+ psin 0) })dpdG \ 
r 77-/2 fx re r-x/cose r 2* r-y /sine \ 

+ + - (U + pcose) 

Jo Jo Jjt/2 Jo Jo Jo / 

exp {— p(A..r ; cos + k u . sin 0) })DpdO 



= k Xj X yj e\p{-x\ Xj -yk y . } [N u (x) + N 2J (x) + N v (x) ] . 
Now let 



(32) 



r "- x i, + x i, 



COS <Pj = kxJfj 



* = e - *j 

sin <t>, = ky.lrj. 



Then 



j(x)= ^ J " U + pcosB) exp {-po eos (0-<I\,)} dpdG 

p/ 2 r/e r nl2 cos ed e 

~* Jo (>< >s (9-0;) + J () rj cos-' (9-*j) 

^ f nl 2 cl,; ^ J^ f w/ 2 ^J COS ( V + <&j ) dV 

cos M' rj J d>j cos 2x P 

i r f*72<u 

sec ^d^-h^ 






x k Xj \ 



f77/2'|) r7T/2<l>. 

4> ; y sec MM^- sin *j y 



sec M'tan *d* 



-+-i J Hog 



^*j / r j~^~ ^- r J 



^yj \ r J — ^i 



+ wr (x «~ x *' ) 



and 



re /* j-/ cose 
2j(*) = I (x + pcos 6) exp {— p(X xj cos0+Xyj sin0)} dpd0 

J W2 J() 



= # 



° 1 -exp fa cos (e-Oj! 



cos 



r jC os (0-*j) 



rf0 



'7T/2 



COS0 1- 



eXP » ^s~0 rjC ° S (e_(|) ^ 



cos 



rj cos (0-$,) + ! 



dQ 



lrr/2 



rj cos 2 (0-Oj) 



211-694 O - 76 - 5 
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CItt r-y/sinQ 

N 3 j(x) = I (x + p cos B) exp { — p (\ x . cos B -f X v . sin B ) }dpdQ 

Je Jo J J 



-.£.-■ 

Jo 



^ 1 - exp -7—^ o cos (B - <I), 



sin B 



j: 



cos 6 1 



rj cos (B - <£>, ) 

y 



rfe 



exp 



in B 



r,- cos (B — <&j) 



y 



sin B 



r/cos (0-4>j)+l 



rfB 



rj cos 2 (B — Oj) 

Entirely analogous calculations give 

/ r 77-/2 fx rG r-xjcosQ ri-n r-y/sinQ 

D 'M-**ML I + L Jo + J, I 

(exp { — k. r . (x 4- p cos B ) — X (/ . (y 4- p sin B ) } ) dpdQ J 



X Xj X yj exp { - sAxj - yX tfj } I -log 



A. j*. / Tj i A. j- 



+ 



J 77/2 



e ! ~ exp — — - r, cos (B - % 



cos 6 



rj cos (B - <&, ) 



^!l ; \ r j \»/ 



dB 



j; 



- 1 " exp ^ shTe ° cos (e_ *' ) 



/•j cos (9 — <f>; 



rf0) 



In order to evaluate \\i (x) use the notation in (31). where now 

zj(x) = pjk XJ k Wj \ | [(x - {)■>+ (y-vr] " 2 exp {- (X,/ + \„jij) }<*£*, 

•>o Jo 

= /3jX Xj .X yj .exp {- «X, r yxJjT « f * + P f =* + ('" j^e 

IaJ o J o J - J o J tf J o 

(p-exp{— p(X J - j cosB-h X /yj sinB)})dpdB] 
=^jX Xj X yj .exp {- zX.*.- yX //} } [zy(*) 4 z 2 j(x) 4- zsjO*)]. 



Now we have 



Jo i A ^j e 



2dO 



osB+X Wj sinB) : 



2 i^ 



de 



: + — . +— Jog 



rj \Xj,j X x j/ rj 



»co8 3 (e-*j) r 

X,, /O + X 



TT 



*j 



sec'^M^ 



; / j *j 



[_Xyj Vj Xyj 



(33) 



(34) 
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Z-ij(x) 



Z3j(x) 



exp rj cos (8- 

1 cost) 



<p;)} 



cosB 



rj cos (B-Oj) + 1) 2 + 1 \dO 



- r 

J e 



exp \-r~7.rj cos (() — *;) 



r?cos 3 (0-*j) 

-y 



cos6 



o cos (e-*,-) + 1 2 + i 



do 



rj cosHQ- <&j) 



5. The Bivariate Symmetric Exponential Density Function 

In this case all the double integrals appearing in (1) and (2) can be reduced to single integrals. 
Using the notation in (14) we now have 



ZC IT- 



Nj(x) = (l/4)A,./, /; §[(x - £)'-' + (y - vY 1 ] " ,/2 exp {- \ XJ \$ - 9 Xj \- k yj \ V - yj \}d^d V 



x J — x 



Dj(x)= (ll4)\ Xj k v \'[\ [(*-£)* + (y-rt) 1 ] -w exvU-\ xi \(-6 x .\-\ x .\y)-6 vj \}dtdT). 
Consider only the first component of Nj{x ) which we again denote by Nj(x). We have 

N J (x)=(m\ Xj ky j n e ' J j* £\x-i\ ' exp {-k Xj (6xj-{)-kyj{6 ¥j -Ti)}d€di l 

+ J J f l*-£l 'exp {-\ Xj (8 x .-t)-\ yj (r,-6 yj )}d£dr) 
+ j ' 'j*€\*-€\ ' ex P {-^/f-«xP-^Wir;-' ri)}d£dit 

t || *l*-£| 'exp {- X*/f — fl*;) - X.jCti - Vj ) }<I^Y| ) 

= (ll4 f )k X jkyj(Nij{x)+N 2 j(x)+N 3 j{x)+N^(x)). 

Now let g — x = p cos 0, y)-y=p sin 0, R 2 j= (d Xj — x) 2 + (d y . — y) 2 , cos 0j = {6 x -x)IRj, sin 
@j = (Oh ~ y)IR i- The limits of integration are determined differently depending upon whether 
(cos Oj, sin 0;) is in the first, second, third, or fourth quadrant. We illustrate the procedure by 
calculating N\j(x) assuming that (cos 0j, sin Oj) is in the first quadrant. We then have 



^^)=exp{-x^ J .-,)-x SJ (^-y)}-{/_j;^ + /7 o Js +/_;/; 

((i + pcos 0) exp {p(k. r cos 0-h \ y . sin 0)})dpd6 
Now let ry= \-,..+ \* cos <l> ? = X r ./rj , sin <&,•= X,* ./r, to obtain 
N iJ (x) = e\p{-\ Xj (d XJ -x)-\y j (e U j-y)} 

i r()j fRjCosdj/cose r n r/jjsinflj/sinfl r-77/2 fx 
I J-tt/2 Jo J«j Jo J -77 Jo 

((% + p cos 0) exp {prj cos (0 — ^j)}) dpdd. 
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The integration with respect to p can be performed explicitly in each term to give 
N lj (x) = exp{-k Xj (d Xj -x)-ky j (0 y -y)} 




\ rj cos (0-<S>j) J 



cos 6 



cos 6j 



Rj cos Oj 



^OS U / I <JU!3 Uj \ r f\j COS Vj }\ 

+ r ]c „ S H0 -*,) ( l -{- Rj ^6 nCOS W"*i) + 1 ) ex P{ ^o7^-^ C ° S < fl -^)j) 



r| cos 2 (0 — O,) 




/ exp Kj — : — — rj 
x\ _ [ sin J 

\ rj cos ( 6 — <£>j ) 



cos (0 — <&j) 



.^-^^ (^Mi^T-Ocos (fl-frj) + l) expf flj^f O-cos (g-ft,) 



sin 



sin 









dd 



de 



(35) 



J J 



The last two terms have been evaluated by transforming the range of integration to 



0.J 



and 



then proceeding as in (33) and (34). Similar expressions are obtained for N 2 j(x), N 3 j(x), N 4j (x) 
and for the analogous terms in Dj(x). The objective function \ff(x) in (1) is evaluated analogously 
except that the integrand with respect to p is of the form p 2 exp (Ap), with A independent of p. 



6. Numerical Results and Discussion 



For (f)j(xj) given by (A), (B), and (C) we have performed a large number of numerical compu- 
tations to find the minimizing P, and \p. In [6], we presented an example illustrating the locally 
linear convergence to P for the case (A), bivariate normal density functions. We do the same here 
for the case (B), exponential density functions, and the case (C), symmetric exponential density 
functions. We also show how the use of Steffensen's iteration [4], accelerates the convergence. 
In order to emphasize the contrast between this probabilistic problem and the corresponding 
deterministic one, we indicate those cases in which x $ co {jlk} £ =1 with an asterisk. Comparisons 
are also made between i//, the bound for \\i from (4), and i// BOUND from (5). 

Problem No. 1 below exhibits a typical set of data used in one of a series of problems studied 
numerically for case (B), the bivariate exponential density function. In this problem A^=20, and the 
4>j(xj) for j= 1,2, . . ., N were taken as bivariate exponential density functions whose values of 
k X j and X yj were generated randomly. The weights f3j were also generated randomly. The iteration 



was terminated when I x n +\ — x n 



io- 



Yn+l — yn 



10 4 . xo was chosen to be 



.7=1 
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where a, 1 = (A !, A,,- 1 ) . M in table 3 is the number of iterations required to achieve the desired 

■i v .1 i Vj > 

accuracy. The computer used for all calculations was the CDC Cyber 72. 



IVoUrm No. 1 



Xx, : 


1.36 


3.56 


13.59 


10.33 


6.30 


5.37 


15.76 


16.43 


4.78 


8.94 




6.68 


17.18 


7.39 


15.90 


0.31 


7.03 


19.88 


12.86 


14.85 


10.94 


K) : 


7.97 


16.21 


8.78 


8.82 


15.01 


2.50 


0.18 


10.35 


18.61 


8.07 




18.00 


7.36 


7.87 


2.90 


9.40 


11.18 


10.76 


6.54 


13.61 


12.66 


to- 


9.06 


8.76 


9.85 


7.94 


2.82 


8.83 


2.39 


2.47 


9.41 


5.70 




7.36 


9.00 


0.94 


0.50 


1.27 


0.13 


6.49 


5.72 


4.28 


9.22 



The method of numerical solution employed is to use (2). This is discussed in more detail in [6]. 
Here we simply wish to illustrate the same kind of linear local convergence for the bivariate ex- 
ponential density function as was shown in |6| for the bivariate normal density function. In table 2, 
we show this local convergence. 

Table 2. Local convergence— bivariate exponential density function 



n 


X n 


y« 


x„ — a 


y n -p 


|*„-a| 


l*.+.-al 


1 


0.326081 


0.144027 


-0.071494 


0.049108 


0.08673 


0.5442 


2 


.354470 


.114142 


-.043105 


.019223 


.04720 


.5983 


3 


.370288 


.102158 


- .027292 


.007239 


.02824 


.6317 


4 


.379962 


.097750 


-.017613 


.00283 1 


.01784 


.6373 


5 


.386272 


.096122 


-.011303 


.001203 


.01 136 


.6356 


6 


.390370 


.095191 


- .007205 


.000572 


.00722 


.6191 


7 


.393015 


.095225 


- .004560 


.000306 


.00147 


.6398 


8 


.394707 


.095100 


- .002868 


.000181 


.00286 


.6049 


9 


.395790 


.095028 


-.001785 


.000109 


.00173 


.6300 


10 


.396482 


.094985 


- .001093 


.000066 


.00109 


.5963 


II 


.396925 


.094958 


- .000650 


.000039 


.00065 


.5692 


12 


.397207 


.094941 


- .000368 


.000022 


.00037 


.5135 


L3 


.397387 


.091930 


-.000188 


.00001 1 


.00019 


.3684 


14 


.397502 


.094923 


- .000073 


.000004 


.00007 





a= (0.397575, 0.094919). 

In table 3 we indicate for 40 randomly generated problems the values of i//, the bound on i// 
from (4), ijj BOUND from (5), and whether or not x is in the convex hull of the means. An asterisk 
next to the value of i// indicates that x is not in the convex hull of the means of distributions. The 
closeness of the bound is probably strongly dependent upon the nature of the data. For these ran- 
domly generated problems, the bounds are not very close. 

TABLE 3. Hounds on objective function — exponential density 



N 


M 


* 


Bound for i// 
from (4) 


BOUND 

from (5) 


30 


16 


29.114 


173.880 


242.628 


30 


12 


12.079 


100.406 


140.618 


30 


14 


5.284* 


55.537 


75.538 


30 


14 


25.792 


202.893 


284.551 


30 


6 


12.656 


77.019 


106.534 


30 


13 


21.548 


347.880 


490.120 


30 


12 


17.742 


118.154 


165.042 


30 


10 


9.856 


81.457 


112.372 


30 


10 


7.350 


61.793 


83.795 
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Table 3. Bounds on objective function — exponential 
density— Continued 



n 


M 


* 


Bound for i// 
from (4) 


ifj BOUND 

from (5) 


30 


10 


8.216 


79.363 


110.170 


20 


14 


5.927 


53.949 


74.012 


20 


14 


26.725 


151.732 


212.956 


20 


14 


8.522 


68.464 


95.960 


20 


11 


3.125 


39.150 


53.083 


20 


20 


23.835 


168.143 


236.237 


20 


10 


4.149 


51.055 


70.417 


20 


8 


9.196* 


54.956 


75.920 


20 


14 


22.643 


344.637 


486.403 


20 


9 


2.400 


25.562 


34.767 


20 


6 


16.436 


105.702 


147.579 


10 


10 


1.121* 


12.074 


15.948 


10 


11 


3.318* 


21.399 


29.313 


10 


9 


0.989* 


12.832 


16.929 


10 


10 


5.497* 


47.132 


65.906 


10 


9 


1.624 


20.033 


27.293 


10 


14 


1.382* 


13.027 


17.314 


10 


10 


4.314 


28.651 


39.116 


10 


14 


5.008 


43.205 


59.725 


10 


15 


1.140* 


10.676 


14.489 


10 


13 


1.773 


24.635 


34.190 


5 


15 


2.271* 


14.688 


20.107 


5 


9 


1.670* 


26.527 


36.952 


5 


7 


0.240* 


3.334 


4.015 


5 


13 


21.036 


106.369 


150.094 


5 


18 


2.140 


13.529 


18.449 


5 


8 


1.075* 


18.091 


25.193 


5 


14 


2.423* 


13.535 


18.639 


5 


28 


4.310 


26.844 


37.449 


5 


17 


1.158* 


12.988 


16.715 


5 


15 


1.438* 


8.857 


12.148 



Problem No. 2 below exhibits a typical set of data for the problems of case (C), the bivariate 
symmetric exponential function. In this problem iV=10, and the (f)j(xj) for j= 1, 2, . . ., TV were 
taken as bivariate symmetric exponential density functions whose values of k, r ., \ y ., 6 Xj , 9 Vj were 
generated randomly. The weights /3j were similarly generated. The iteration was terminated when 
\x n +i — x n \ ^ 10 4 and |y n +i ~Jn\ ^ 10 4 . ^o was chosen to be: 

£ m 

j=i 






where dj= (6x- r Uj ). M in table 5 is the number of iterations required to achieve the desired 
accuracy. 

Problem No. 2 

k x .: 1.36 17.53 8.82 5.37 4.78 18.61 6.68 18.00 2.90 7.03 
j 

X«.: 7.97 13.59 15.89 2.50 16.43 18.82 18.00 7.39 1.01 11.18 

6 Xj : 9.06 4.39 3.15 8.83 5.18 4.47 7.36 3.94 0.15 0.13 

Byy 1.78 9.85 7.51 7.88 2.47 4.04 8.59 0.94 4.70 9.94 

Pj: 8.11 5.17 2.82 0.09 2.39 5.70 3.68 7.95 1.27 5.38 
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Local convergence for Problem No. 2 is exhibited in table 4. In table 5 we indicate the bounds and 
whether or not the solution is in the convex hull of the means in the same fashion as we did for the 
bivariate symmetric exponential density {'unction. It is similar to table l\. For the problems generated 
for the symmetric exponential case we never found one in the convex hull of the means. 
Table 4. Local convergence— bivariate symmetric exponential density function 



n 


X„ 


y» 


x n — a 


y»-fi 


\x n -a\ 


|x/m t-5| 
\x n -a\ 


1 


0.01668185 


0.04066205 


0.00462450 


0.03773900 


0.0380213 


0.3519 


2 


.01158224 


.01629392 


-0.O0O47511 


.01337087 


.0133793 


.3051 


3 


.01210797 


.00700439 


.00005062 


.00408134 


.00408165 


.3864 


4 


.01255962 


.00441807 


.00050227 


.00149502 


.00157714 


.5018 


5 


.01253156 


.00355670 


.00047421 


.00063365 


.000791446 


.5563 


6 


.01238560 


.00321650 


.00032825 


.00029345 


.000440296 


.5251 


7 


.01224484 


.00305836 


.00018749 


.00013531 


.000231217 


.4034 


8 


.01213578 


.00297353 


.00007843 


.00005048 


.0000932711 





a = (0.0 1 205735. 0.00292305). 

Table 5. Hounds on objective function — symmetric exponential density 













N 


M 


* 


Hound for i// 
from (4) 


\\f Hound 
from (5) 


30 


13 


38.195* 


185.942 


533.053 


30 


19 


36.889* 


609.044 


699.347 


30 


II 


31.590* 


623.526 


696.993 


30 


14 


15.391* 


513.657 


569.532 


30 


13 


26.269* 


850.293 


978.537 


30 


11 


3.362* 


527.453 


596.923 


30 


13 


15.910* 


754.403 


823.259 


30 


28 


54.418* 


553.564 


611.917 


30 


10 


56.257* 


723.158 


801.155 


30 


12 


7.493* 


182.261 


544.145 


20 


16 


0.02035* 


324.837 


349.831 


20 


17 


119.638* 


503.372 


565.895 


20 


10 


5.548* 


312.673 


350.797 


20 


12 


10.965* 


369.864 


420.667 


20 


16 


17.922* 


469.078 


527.044 


20 


25 


199.676* 


682.171 


762.248 


20 


13 


72.808* 


553.628 


631.301 


20 


21 


2.169* 


215.444 


241.296 


20 


9 


4.340* 


452.851 


483.389 


20 


13 


0.45620* 


342.029 


376.407 


10 


8 


2.407* 


178.300 


193.276 


10 


17 


29.185* 


308.299 


366.341 


10 


15 


4.695* 


107.535 


117.199 


10 


9 


0.47852* 


177.145 


187.936 


10 


12 


5.742* 


237.441 


254.580 


10 


19 


25.345* 


205.296 


240.394 


10 


12 


9.546* 


186.705 


203.622 


10 


12 


1.429* 


140.736 


150.786 


10 


20 


2.901* 


155.144 


1 78.605 


10 


5 


0.77443* 


286.310 


307.41 1 


5 


12 


9.236* 


93.932 


124.836 


5 


25 


0.00133* 


49.655 


55.333 


5 


4 


1.951* 


85.311 


94.246 


5 


19 


14.127* 


216.730 


273.959 


5 


8 


0.00000* 


88.316 


93.744 


5 


8 


2.465* 


60.637 


67.276 


5 


17 


3.493* 


100.387 


110.473 


5 


13 


1.301* 


80.547 


85.485 


5 


23 


30.839* 


69.109 


79.990 


5 


13 


0.00000* 


31.311 


35.560 



71 



6.1. Steffensen's Iteration 

A well known method for accelerating the convergence of linearly convergent iterative schemes 
is Steffensen's iteration [4]. In the case when the unknown vector is two dimensional, Steffensen's 
iteration, applied to the scheme in (2), is described as follows (see [4], p. 115 for a derivation): 



Choose a vector x (0) = (x 



(0) 



r 



(0) 



), and construct the sequence of vectors {x ik) } by the rule: 



for each A = 0, 1, 2, . . . , set xo = x {k \ and calculate X\ , x 2 , X3 from 



Xn - 



H(x„) 



0, 1, 2. 



Let the matrix X„ be defined by 



Xn = 



Xn X n + 1 



\Yn Jn + U 

and let A denote the forward difference operator. The next iterate is given by 

*<* + 1 > = *q-AXo(A 2 *o)- 1 A* . 

In order to test the utility of Steffensen's iteration in reducing the amount of calculation to 
converge to any given degree of tolerance, ten, 20 destination problems and ten, 10 destination 
problems were run on the CDC Cyber 72 both with and without the use of Steffensen's iteration 



Table 6. Effect of Steffensen's iteration on convergence 



Number of 


Number of 


Number of iterations 


destinations 


iterations 


(Steffensen) 


20 


14 


7 


20 


14 


13 


20 


16 


6 


20 


16 


13 


20 


17 


15 


20 


14 


6 


20 


12 


5 


20 


15 


5 


20 


15 


9 


20 


13 


6 


20 


14 


6 


10 


15 


12 


10 


15 


7 


10 


16 


6 


10 


16 


8 


10 


18 


9 


10 


12 


10 


10 


17 


10 


10 


13 


8 


10 


25 


13 


10 


16 


7 



technique. The number of iterations given in table 6 is the total number of iterations, i.e., the num- 
ber reported using Steffensen's iteration includes the three iterations required to use the method. 
The problems all used a bivariate normal density function. The average ratio of the number of itera- 
tions with the Steffensen method to the number of iterations without the method for the 20 problems 
of table 6 was 0.571. Hence, it can be readily seen that it is quite advantageous to use the Steffensen 
iteration method when performing these calculations. 



We thank the reviewer for his helpful suggestions and for his careful verification of many 
of the formulas presented in this paper. 
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